rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-figure_a14a
#-figure_a14b
#-figure_a14c
#-figure_a14d

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library("ggrepel")
library("lemon")
library("sandwich")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("pgirmess")
library("geojsonsf")

#Step1. Turning some variables to numeric
data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)
data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)

data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
data_cov_mob$bridging_tot_pop_dummy<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)


data_cov_mob_week<-subset(data_cov_mob, week=="01")
data_cov_mob_week<-subset(data_cov_mob_week, select = c("GEN",
                                                        "bridging_tot_pop"))

data<-data_cov_mob
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)
data$post_treat_bridging_tot_pop_dummy<-data$post_treat*data$bridging_tot_pop_dummy

###########
#All clubs#
###########


#Create dummies

library('fastDummies')
# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')

library(dplyr)
new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))

list_dummies<-names(new)
#I am removing the counties which get dropped in Stata

#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]
#Re-enable this later if necessary
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_09"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_52"]


list_land_dummies<-names(new3)
#Re-enable this later if necessary
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]


#Event study
#Bridging Interaction
new2<-dataw %>% dplyr::select(starts_with("week_"))
list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]


############
#change_afd#
############


#Step1 Creating the interaction terms for change_afd_sec_narrow
for (i in list_week_dummies){
  datas[[paste("change_afd_", i, sep="")]]<-datas[[i]]*datas$pct_change_afd
}


#Step2 Selecting the right independent vars
inter_terms_change_afd<-datas %>% dplyr::select(starts_with("change_afd_week"))
inter_terms_change_afd <-names(inter_terms_change_afd)
inter_terms_change_afd <-inter_terms_change_afd[inter_terms_change_afd != "change_afd_week_10"]

#Deleting variables which are collinear
list_land_dummies_full2<-list_land_dummies_full[list_land_dummies_full!= "SN_L_02"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_53"]

#Step3 Creating the list of variables
for_event_inter_terms_change_afd<-list.append(inter_terms_change_afd, list_land_dummies_full, list_week_dummies)
for_event_inter_terms_change_afd2<-list.append(inter_terms_change_afd, list_land_dummies_full2, list_week_dummies2)

#Step4 Deleting missing observations
datas$time<-as.numeric(datas$week)
case_change_afd<-subset(datas, select = list.append(for_event_inter_terms_change_afd2, "pct_afd", "mean_mobil", "POINT_X", "POINT_Y",
                                                    "time", "AGS"))
case_change_afd <-case_change_afd[complete.cases(case_change_afd), ]


#Step5 Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(case_change_afd$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Step6 Merging with spatial dataframe
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
case_change_afd_oneweek<-subset(case_change_afd, week_11==1)

#Step7 Merging with shape file
nuts_merged<-left_join(nuts2, case_change_afd_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$pct_afd))

#Step8 Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Step9: Calculating the correlogram
#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$pct_afd, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
dist_cut_change_afd<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

#Step10: change_afd formula
time_trend_change_afd<-as.formula(paste("mean_mobil~",
                                        paste(for_event_inter_terms_change_afd2, collapse= "+")))
dm_change_afd <- dist_mat(case_change_afd, lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time")

#Step11: Run the Conley SE
#Takes 30 sec
df_change_afd<-conleyreg(time_trend_change_afd, data=case_change_afd,
                         dist_cut_change_afd, unit = "AGS", time = "time", 
                         dist_mat=dm_change_afd, lag_cutoff=temp_lag)

df_change_afd<-broom::tidy(df_change_afd)
df_change_afd$term[grepl( "change_afd_week_" , df_change_afd$term)]
df_change_afd<-add_row(df_change_afd,term="change_afd_week_10",estimate=0)
df_change_afd<-subset(df_change_afd, grepl( "change_afd_week_" , df_change_afd$term))
df_change_afd<-df_change_afd[order(df_change_afd$term),]


#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_change_afd$observation<-NA
df_change_afd$observation <- x

df_change_afd$observation<-as.character(df_change_afd$observation)
selection<-df_change_afd$observation[!stringr::str_detect(df_change_afd$observation, '-')]
df_change_afd$observation[!stringr::str_detect(df_change_afd$observation, '-')]<-paste("+", selection, sep = "")
df_change_afd$t<-paste("Pct. AfD x \n", "(t", df_change_afd$observation, ")", sep = "")
df_change_afd$t[df_change_afd$t=="Pct. AfD x \n(t+0)"]<-"Pct. AfD x t"

df_change_afd$week<-NA
df_change_afd$week<-readr::parse_number(df_change_afd$term)
df_change_afd$model<-"Pct. AfD"


#################
#All Event Study#
#################
irislabs1<-seq(1,53,6)
irislabs2<-unique(df_change_afd[as.numeric(df_change_afd$week) %in% irislabs1,]$t)

coef_change_afd<-ggplot(df_change_afd, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_change_afd, file = "paper/graphs/figure_a14a.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



###########################
#Zoom Event Study Bridging#
###########################

df_change_afd_small<-subset(df_change_afd, week>=(5) & week<=(16))

irislabs1<-seq(1,53,1)
irislabs2<-unique(df_change_afd[as.numeric(df_change_afd$week) %in% irislabs1,]$t)

coef_zoom_change_afd<-ggplot(df_change_afd_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_change_afd, file = "paper/graphs/figure_a14c.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



##############
#AFD as Dummy#
##############

mean_pct_afd<-summary(datas$pct_change_afd)[4]
datas$pct_change_afd_dummy<-ifelse(datas$pct_afd>mean_pct_afd, 1, 0)

for (i in list_week_dummies){
  datas[[paste("change_afd_dummy_", i, sep="")]]<-datas[[i]]*datas$pct_change_afd_dummy
}


#Step2 Selecting the right independent vars
inter_terms_change_afd_dummy<-datas %>% dplyr::select(starts_with("change_afd_dummy_week"))
inter_terms_change_afd_dummy <-names(inter_terms_change_afd_dummy)
inter_terms_change_afd_dummy <-inter_terms_change_afd_dummy[inter_terms_change_afd_dummy != "change_afd_dummy_week_10"]

#Deleting variables which are collinear
list_land_dummies_full2<-list_land_dummies_full[list_land_dummies_full!= "SN_L_02"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_53"]

#Step3 Creating the list of variables
for_event_inter_terms_change_afd_dummy<-list.append(inter_terms_change_afd_dummy, list_land_dummies_full, list_week_dummies)
for_event_inter_terms_change_afd_dummy2<-list.append(inter_terms_change_afd_dummy, list_land_dummies_full2, list_week_dummies2)

#Step4 Deleting missing observations
datas$time<-as.numeric(datas$week)
case_change_afd_dummy<-subset(datas, select = list.append(for_event_inter_terms_change_afd_dummy2, "pct_afd", "mean_mobil", "POINT_X", "POINT_Y",
                                                          "time", "AGS"))
case_change_afd_dummy <-case_change_afd_dummy[complete.cases(case_change_afd_dummy), ]


#Step5 Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(case_change_afd_dummy$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Step6 Merging with spatial dataframe
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
case_change_afd_dummy_oneweek<-subset(case_change_afd_dummy, week_11==1)

#Step7 Merging with shape file
nuts_merged<-left_join(nuts2, case_change_afd_dummy_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$pct_afd))

#Step8 Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Step9: Calculating the correlogram
#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$pct_afd, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
dist_cut_change_afd_dummy<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

#Step10: change_afd_dummy formula
time_trend_change_afd_dummy<-as.formula(paste("mean_mobil~",
                                              paste(for_event_inter_terms_change_afd_dummy2, collapse= "+")))
dm_change_afd_dummy <- dist_mat(case_change_afd_dummy, lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time")

#Step11: Run the Conley SE
#Takes 30 sec
df_change_afd_dummy<-conleyreg(time_trend_change_afd_dummy, data=case_change_afd_dummy,
                               dist_cut_change_afd_dummy, unit = "AGS", time = "time", 
                               dist_mat=dm_change_afd_dummy, lag_cutoff=temp_lag)

df_change_afd_dummy<-broom::tidy(df_change_afd_dummy)
df_change_afd_dummy$term[grepl( "change_afd_dummy_week_" , df_change_afd_dummy$term)]
df_change_afd_dummy<-add_row(df_change_afd_dummy,term="change_afd_dummy_week_10",estimate=0)
df_change_afd_dummy<-subset(df_change_afd_dummy, grepl( "change_afd_dummy_week_" , df_change_afd_dummy$term))
df_change_afd_dummy<-df_change_afd_dummy[order(df_change_afd_dummy$term),]


#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_change_afd_dummy$observation<-NA
df_change_afd_dummy$observation <- x

df_change_afd_dummy$observation<-as.character(df_change_afd_dummy$observation)
selection<-df_change_afd_dummy$observation[!stringr::str_detect(df_change_afd_dummy$observation, '-')]
df_change_afd_dummy$observation[!stringr::str_detect(df_change_afd_dummy$observation, '-')]<-paste("+", selection, sep = "")
df_change_afd_dummy$t<-paste("Pct. AfD x \n", "(t", df_change_afd_dummy$observation, ")", sep = "")
df_change_afd_dummy$t[df_change_afd_dummy$t=="Pct. AfD x \n(t+0)"]<-"Pct. AfD x t"

df_change_afd_dummy$week<-NA
df_change_afd_dummy$week<-readr::parse_number(df_change_afd_dummy$term)
df_change_afd_dummy$model<-"Pct. AfD"



#################
#All Event Study#
#################
irislabs1<-seq(1,53,6)
irislabs2<-unique(df_change_afd_dummy[as.numeric(df_change_afd_dummy$week) %in% irislabs1,]$t)

coef_change_afd_dummy<-ggplot(df_change_afd_dummy, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_change_afd_dummy, file = "paper/graphs/figure_a14b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)





###########################
#Zoom Event Study Bridging#
###########################

df_change_afd_dummy_small<-subset(df_change_afd_dummy, week>=(5) & week<=(16))

irislabs1<-seq(1,53,1)
irislabs2<-unique(df_change_afd_dummy[as.numeric(df_change_afd_dummy$week) %in% irislabs1,]$t)

coef_zoom_change_afd<-ggplot(df_change_afd_dummy_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_change_afd, file = "paper/graphs/figure_a14d.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



